library(xlsx)
library(plyr)
library(dplyr)
library(tidyr)
library(ARTool)
library(ggplot2)
library(ggsignif)
library(ggpubr)
library(gridExtra)

setwd("H://OneDrive//研究//3-小实验TU//5-论文//TU对本外竞争的影响//数据与分析")
                                                                           #设置数据读取与输出文件夹

##################
###  数据输入  ###
##################
rm(list=ls())
grow <- read.xlsx("20220712.xlsx",sheetName = "grow", encoding = "UTF-8")  #生长数据
grow$Treat <- as.factor(grow$Treat)
grow$W <- as.factor(grow$W)
grow$U <- as.factor(grow$U)
grow$ID <- as.factor(grow$ID)
grow$PT <- as.factor(grow$PT)

leaf <- read.xlsx("20220712.xlsx",sheetName = "leaf", encoding = "UTF-8")  #叶片数据
leaf$W <- as.factor(leaf$W)
leaf$U <- as.factor(leaf$U)
leaf$ID <- as.factor(leaf$ID)
leaf$PT <- as.factor(leaf$PT)

ps <- read.xlsx("20220712.xlsx",sheetName = "ps", encoding = "UTF-8")     #光合数据
ps$W <- as.factor(ps$W)
ps$U <- as.factor(ps$U)
ps$ID <- as.factor(ps$ID)
ps$PT <- as.factor(ps$PT)

cnp <- read.xlsx("20220712.xlsx",sheetName = "cnp", encoding = "UTF-8")  #元素数据
cnp$W <- as.factor(cnp$W)
cnp$U <- as.factor(cnp$U)
cnp$ID <- as.factor(cnp$ID)
cnp$PT <- as.factor(cnp$PT)


##################
###  定义函数  ###
##################
summary.yb <- function(data, group, value, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {
  library(plyr)
  library(dplyr)
  library(ARTool)
  result <- list()
  for (i in c(1:length(value))) {
    result[[i]] <- list()
    length2 <- function (x, na.rm=T) {
      if (na.rm) sum(!is.na(x))
      else       length(x)
    }
    datac <- ddply(data, group, .drop=T,
                   .fun = function(xx, col) {
                     c(N    = length2 (xx[[col]], na.rm=T),
                       mean = mean    (xx[[col]], na.rm=T),
                       sd   = sd      (xx[[col]], na.rm=T),
                       min  = min     (xx[[col]], na.rm=T),
                       max  = max     (xx[[col]], na.rm=T)
                     )
                   },
                   value[i]
    )
    datac$se <- datac$sd / sqrt(datac$N) 
    datac$down <- datac$mean - datac$se
    datac$up <- datac$mean + datac$se
    result[[i]][[1]] <- datac
    
    m <- art(eval(parse(text = paste(paste(value[i],"~"), paste(group, collapse = "*")))), data=data)  
                                                                     #根据需要进行多因素分析
    result[[i]][[2]] <- anova(m)                                     #查看交互分析结果
    result[[i]][[3]] <- art.con(m, eval(parse(text = paste(paste("~ "), paste(group, collapse = "*")))))%>%
      summary() %>%
      mutate(sig = symnum(p.value, corr = FALSE, na = FALSE,
                          cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                          symbols = c("***", "**", "*", ".", " ")
                          )
             )                                                       #查看各单元之间交互作用的具体结果
  }
  return(result)
}   #输入: data-数据框，group-自变量(向量), value-因变量(可多个, 向量);
    #输出: 列表, 一级-因变量(输入顺序), 二级1-统计描述, 二级2-非参数方差分析F, 二级3-非参数多重比较

change.yb <- function(data, group, treat, ck, value, list) {
  library(plyr)
  level <- unique(data[,treat])
  for (i in c(1:length(value))) {
    mean_ck <- ddply(data[data[,treat] == ck,], group, mutate, ck = mean(eval(parse(text=value[i])), na.rm=T))
    dataa <- ddply(data, group, mutate, vra = (eval(parse(text=value[i]))))
    vra <- data.frame()
    for (j in c(1:length(level))) {
      vra <- rbind(vra, 
                   data.frame(Treat = level[j], 
                              vra = (dataa[dataa[,treat] == level[j], "vra"] - mean_ck[,"ck"])/mean_ck[,"ck"]))
    }
    data_vra <- cbind(dataa[,setdiff(group, treat)], vra)
    datab <- ddply(data_vra, group, summarize, N = length(vra), vra_mean = round(mean(vra),4), vra_sd = round(sd(vra),4))
    datab$vra_se <- datab$vra_sd / sqrt(datab$N) 
    datab$vra_down <- datab$vra_mean - datab$vra_se
    datab$vra_up <- datab$vra_mean + datab$vra_se
    list[[i]][[4]] <- datab
  }
  return(list)
}   #输入: data-数据框, group-自变量(向量), treat-相对值计算分组(单一字符串), ck-相对值计算的目标值, value-因变量(可多个, 向量), list-原有结果列表;
    #输出: list[[计算目标值]][[4]]

transpose.yb <- function(list, name) {
  New <- list()
  for (i in c(1:length(list))) {
    for (j in c(1:length(list[[i]]))) {
      dataf <- cbind(data.frame(Index = rep(name[i], length(list[[i]][[j]][,1]))), list[[i]][[j]])
      if (i == 1) {
        New[[j]] <- dataf
      } else {
        New[[j]] <- rbind(New[[j]], dataf)
      }
    }
  }
  return(New)
}     #将多重列表转化为一层列表，用于输出excel

fig.yb <- function(list, fig_name) {
  library(dplyr)
  library(ggplot2)
  library(ggpubr)
  library(gridExtra)
  library(do)
  
  fig_index <- list()
  
  for (i in c(1:length(list))) {
    actual_data <- list[[i]][[1]]
    actual_data$Treat <- paste(actual_data$W, actual_data$U, sep = "-")
    actual_data$Treat[actual_data$Treat == "0-0"] <- "CK"
    actual_data$Treat[actual_data$Treat == "1-0"] <- "T"
    actual_data$Treat[actual_data$Treat == "0-1"] <- "U"
    actual_data$Treat[actual_data$Treat == "1-1"] <- "T*U"
    actual_data$Treat <- factor(actual_data$Treat, levels = c("CK", "T", "U", "T*U"))
    actual_data$ID <- factor(actual_data$ID, levels = c("Mono", "Mix"))
    
    diff_data <- separate(list[[i]][[3]], contrast, 
                          into = c("T1","U1", "ID1", "PT1", "T2","U2", "ID2", "PT2"), remove = FALSE)
    diff_data <- diff_data[(diff_data$T1 == diff_data$T2 & 
                              diff_data$U1 == diff_data$U2 & 
                              diff_data$p.value < .05) & 
                             (diff_data$ID1 == diff_data$ID2 | 
                                diff_data$PT1 == diff_data$PT2), ]
    diff_data$Treat <- paste(diff_data$T1, diff_data$U1, sep = "-")
    diff_data$Treat[diff_data$Treat == "0-0"] <- "CK"
    diff_data$Treat[diff_data$Treat == "1-0"] <- "T"
    diff_data$Treat[diff_data$Treat == "0-1"] <- "U"
    diff_data$Treat[diff_data$Treat == "1-1"] <- "T*U"
    diff_data$Treat <- factor(diff_data$Treat, levels = c("CK", "T", "U", "T*U"))
    if (length(diff_data$contrast) > 0) {
      for (j in c(1:length(diff_data$contrast))) {
        if (diff_data$ID1[j] == diff_data$ID2[j]) {
          diff_data$y[j] <- max(actual_data$up[actual_data$Treat == diff_data$Treat[j] &
                                                 actual_data$ID == diff_data$ID1[j] &
                                                 actual_data$PT == diff_data$PT1[j]],
                                actual_data$up[actual_data$Treat == diff_data$Treat[j] &
                                                 actual_data$ID == diff_data$ID2[j] &
                                                 actual_data$PT == diff_data$PT2[j]]) +
            .06 * (max(actual_data$up) - min(actual_data$down))
        } else {
          diff_data$y[j] <- max(actual_data$up[actual_data$Treat == diff_data$Treat[j] &
                                                 actual_data$ID == diff_data$ID1[j] &
                                                 actual_data$PT == diff_data$PT1[j]],
                                actual_data$up[actual_data$Treat == diff_data$Treat[j] &
                                                 actual_data$ID == diff_data$ID2[j] &
                                                 actual_data$PT == diff_data$PT2[j]]) +
            .12 * (max(actual_data$up) - min(actual_data$down))
        }
      }
      diff_data$ID1[diff_data$ID1 == "Mono"] <- 1
      diff_data$ID2[diff_data$ID2 == "Mono"] <- 1
      diff_data$ID1[diff_data$ID1 == "Mix"] <- 2
      diff_data$ID2[diff_data$ID2 == "Mix"] <- 2
      diff_data$PT1[diff_data$PT1 == "N"] <- -.125
      diff_data$PT2[diff_data$PT2 == "N"] <- -.125
      diff_data$PT1[diff_data$PT1 == "I"] <- .125
      diff_data$PT2[diff_data$PT2 == "I"] <- .125
      diff_data$ID1 <- as.numeric(diff_data$ID1)
      diff_data$ID2 <- as.numeric(diff_data$ID2)
      diff_data$PT1 <- as.numeric(diff_data$PT1)
      diff_data$PT2 <- as.numeric(diff_data$PT2)
      diff_data$xmin <- pmin(diff_data$ID1 + diff_data$PT1, diff_data$ID2 + diff_data$PT2)
      diff_data$xmax <- pmax(diff_data$ID1 + diff_data$PT1, diff_data$ID2 + diff_data$PT2)
    }
    
    F_data <- as.data.frame(list[[i]][[2]])[as.data.frame(list[[i]][[2]])$`Pr(>F)` < 0.05 ,c("F value", "Pr(>F)")]
    F_data$`F value` <- round(F_data$`F value`, 2)
    F_data$`Pr(>F)` <- round(F_data$`Pr(>F)`, 4)
    F_data$`Pr(>F)`[F_data$`Pr(>F)` < 0.0001] <- "< 0.0001"
    rownames(F_data) <- Replace(rownames(F_data), ":" ," * ")
    rownames(F_data) <- Replace(rownames(F_data), "W" ,"T")
    rownames(F_data) <- Replace(rownames(F_data), "PT" ,"PS")
    rownames(F_data) <- Replace(rownames(F_data), "ID" ,"PT")
    
    vra_data <- list[[i]][[4]]
    vra_data$Treat <- as.character(vra_data$Treat)
    vra_data$Treat[vra_data$Treat == "T0"] <- "CK"
    vra_data$Treat[vra_data$Treat == "T1"] <- "T"
    vra_data$Treat[vra_data$Treat == "U0"] <- "U"
    vra_data$Treat[vra_data$Treat == "U1"] <- "T*U"
    vra_data$Treat <- factor(vra_data$Treat, levels = c("CK", "T", "U", "T*U"))
    vra_data$ID <- factor(vra_data$ID, levels = c("Mono", "Mix"))
    
    for (j in c(1:length(unique(actual_data$Treat)))) {
      if (j == 1) {
        fig_treat <- list()
        vra_fig <- NULL
        
      } else {
        vra_fig <- ggplot(vra_data[vra_data$Treat == c("CK", "T", "U", "T*U")[j],])+
          geom_col(aes(x = ID, y = vra_mean, group = PT, fill = PT), 
                   position = position_dodge(width = -1), width = .5)+
          geom_errorbar(aes(x = ID, ymin = vra_down, ymax = vra_up, group = PT, color = PT),
                        position = position_dodge(width = -1), width = .3)+
          geom_hline(yintercept = 0)+
          ylim(min(vra_data$vra_down), max(vra_data$vra_up))+
          theme_bw()+
          labs(x = NULL, y = NULL)+
          theme(legend.position = "none")
      }
      
      actual_fig <- ggplot(actual_data[actual_data$Treat == c("CK", "T", "U", "T*U")[j],])+
        geom_line(aes(x = ID, y = mean, group = PT, color = PT), linetype =2, position = position_dodge(width = -0.5))+
        geom_pointrange(aes(x = ID, y = mean, ymin = down, ymax = up, color = PT), position = position_dodge(width = -0.5))+
        ylim(min(actual_data$down), max(actual_data$up)*1.12)+
        theme_bw()+
        labs(x = NULL, y = NULL)+
        theme(legend.position = "none")
      if (length(diff_data[diff_data$Treat == c("CK", "T", "U", "T*U")[j],1]) != 0) {
        actual_fig <- actual_fig +
          geom_signif(
            data = diff_data[diff_data$Treat == c("CK", "T", "U", "T*U")[j],],
            aes(xmin = xmin, xmax = xmax, annotations = "", y_position = y), 
            manual = TRUE, vjust = .5
          )+
          geom_text(data = diff_data[diff_data$Treat == c("CK", "T", "U", "T*U")[j],], 
                    aes(x = (xmin + xmax)/2, y = y, label = sig))
      }
      
      fig_treat[[j]] <- ggarrange(actual_fig, vra_fig, nrow = 2, ncol = 1, align = "hv")
      
    }
    fig_index[[i]] <- annotate_figure(
      ggarrange(fig_treat[[1]], fig_treat[[2]], fig_treat[[3]], fig_treat[[4]],
                nrow = 1, ncol = 4, align = "hv"),
      left = text_grob(fig_name[i], size = 12, rot = 90))
  }
  return(fig_index)
}

##################
###  数据分析  ###
##################
#二级数据计算
grow$TB <- grow$AB + grow$UB                                         #总生物量
grow$RS <- grow$UB / grow$AB                                         #根冠比
for (i in c(1:length(grow[,1]))) {
  grow$RY[i] <- grow$TB[i] / mean(grow$TB[grow$Treat == grow$Treat[i] & grow$PT == grow$PT[i] & grow$re == grow$re[i] & grow$ID == "Mono"])
}                                                                    #相对产量
cnp$CN <- cnp$C / cnp$N                                              #碳氮比
cnp$CP <- cnp$C / cnp$P                                              #碳磷比
cnp$NP <- cnp$N / cnp$P                                              #氮磷比
  
#统计描述与非参数方差分析
group <- c("W", "U", "ID", "PT")
result_grow <- summary.yb(grow, group, c("H", "RL", "BS", "TB", "RS", "RY"))
result_grow <- change.yb(grow, c("Treat", "ID", "PT"), "Treat", "T0", c("H", "RL", "BS", "TB", "RS", "RY"), result_grow)
result_leaf <- summary.yb(leaf, group, c("leaf_a", "leaf_m"))
result_leaf <- change.yb(leaf, c("Treat", "ID", "PT"), "Treat", "T0", c("leaf_a", "leaf_m"), result_leaf)
result_ps <- summary.yb(ps, group, c("A", "Ci", "gs", "E", "WUE", "chl", "Fv.Fm"))
result_ps <- change.yb(ps, c("Treat", "ID", "PT"), "Treat", "T0", c("A", "Ci", "gs", "E", "WUE", "chl", "Fv.Fm"), result_ps)
result_cnp <- summary.yb(cnp, group, c("CN", "CP", "NP"))
result_cnp <- change.yb(cnp, c("Treat", "ID", "PT"), "Treat", "T0", c("CN", "CP", "NP"), result_cnp)


##################
###  结果输出  ###
##################
#表格输出
Table <- c(transpose.yb(result_grow, c("H", "RL", "BS", "TB", "RS", "RY")),
           transpose.yb(result_leaf, c("leaf_a", "leaf_m")),
           transpose.yb(result_ps, c("A", "Ci", "gs", "E", "WUE", "chl", "Fv.Fm")),
           transpose.yb(result_cnp, c("CN", "CP", "NP")))
for (i in c(1,5,9,13)) {
  Table[[i]]$Treat <- paste(Table[[i]]$W, Table[[i]]$U, sep = "-")
  Table[[i]]$Treat[Table[[i]]$Treat == "0-0"] <- "CK"
  Table[[i]]$Treat[Table[[i]]$Treat == "1-0"] <- "T"
  Table[[i]]$Treat[Table[[i]]$Treat == "0-1"] <- "U"
  Table[[i]]$Treat[Table[[i]]$Treat == "1-1"] <- "T*U"
  Table[[i]]$Treat <- factor(Table[[i]]$Treat, levels = c("CK", "T", "U", "T*U"))
  Table[[i]]$ID <- factor(Table[[i]]$ID, levels = c("Mono", "Mix"))
}
sheetname <- c("grow_SD", "grow_ANOVA", "grow_MC", "grow_vra", 
               "leaf_SD", "leaf_ANOVA", "leaf_MC", "leaf_vra",
               "ps_SD", "ps_ANOVA", "ps_MC", "ps_vra", 
               "cnp_SD", "cnp_ANOVA", "cnp_MC", "cnp_vra")
for (i in 1:16) {
  if (i == 1) {write.xlsx(Table[[i]], "result.xlsx",  sheetName = sheetname[i])}
  else {write.xlsx(Table[[i]], "result.xlsx",  sheetName = sheetname[i], append = TRUE)}
} 

#绘图
Fig_treat_data <- data.frame("CK" = "Control check (CK)", "T" = "Temperature (T)", "U" = "UV-B (U)", "T × U" = "T × U")
row.names(Fig_treat_data) <- ""
Fig_treat <- tableGrob(Fig_treat_data, theme = ttheme_minimal(core=list(fg_params=list(fontface="bold"))), cols = NULL)
Fig_treat$widths <- unit(c(0.05, rep(0.2375, 4)), "npc")  #绘制处理标注（以表格形式）

Fig_legend_data <- data.frame(
  plant = c("Native", "Invasive"),
  x = c(1, 2),
  y = c(0, 0),
  ymin = c(-1, -1),
  ymax = c(1, 1)
)
Fig_legend <- cowplot::get_legend(
  ggplot(Fig_legend_data)+
    geom_col(aes(x = x, y = y, group = plant, fill = plant))+
    theme_void()+
    scale_fill_discrete(name = "", 
                        breaks = c("Native", "Invasive"), 
                        labels = c("Native", "Invasive"))+
    theme(legend.position="top",
          legend.text = element_text(size = 12, face = "bold"))
)  #绘制图例

Fig_grow <- fig.yb(result_grow, 
                   c("Plant height (cm)", "Root length (cm)", "Basal stem (mm)", 
                     "Total biomass (g)", "Root shoot ratio", "Relative yield"))
Fig_leaf <- fig.yb(result_leaf, 
                   c(expression("Leaf area (" * mm ^ 2 * ")"), 
                     "Leaf fresh weight (mg)")) 
Fig_ps   <- fig.yb(result_ps, 
                   c(expression("Net photosynthetic rate (" * μmol ~ m ^ -2 ~ s ^ -1 * ")"), 
                     expression("Intercellular " * CO[2] * " concentration (ppm)"), 
                     expression("Stomatal conductance (" * μmol ~ m ^ -2 ~ s ^ -1  * ")"), 
                     expression("Transpiration rate (" * μmol ~ m ^ -2 ~ s ^ -1  * ")"), 
                     "Water use efficiency (%)", "Chlorophyll content (SPDA)", "Fv/Fm"))
Fig_cnp  <- fig.yb(result_cnp, c("C:N", "C:P", "N:P"))

Fig1 <- annotate_figure(
  ggarrange(Fig_grow[[1]], Fig_grow[[2]], Fig_grow[[4]], Fig_grow[[5]], Fig_grow[[6]],
            nrow = 5, ncol = 1, align = "hv", common.legend = T, legend = "left"),
  bottom = Fig_legend,
  top = Fig_treat)

Fig2 <- annotate_figure(
  ggarrange(Fig_leaf[[1]], Fig_leaf[[2]],
            nrow = 2, ncol = 1, align = "hv", common.legend = T, legend = "left"),
  bottom = Fig_legend,
  top = Fig_treat)

Fig3 <- annotate_figure(
  ggarrange(Fig_ps[[1]], Fig_ps[[2]], Fig_ps[[3]], Fig_ps[[4]], Fig_ps[[5]], Fig_ps[[6]], Fig_ps[[7]],
            nrow = 7, ncol = 1, align = "hv", common.legend = T, legend = "left"),
  bottom = Fig_legend,
  top = Fig_treat)

Fig4 <- annotate_figure(
  ggarrange(Fig_cnp[[1]], Fig_cnp[[2]], Fig_cnp[[3]],
            nrow = 3, ncol = 1, align = "hv", common.legend = T, legend = "left"),
  bottom = Fig_legend,
  top = Fig_treat)

###存图
ggsave(file = "Fig1.svg", plot = Fig1, width = 10, height = 17.5)
ggsave(file = "Fig2.svg", plot = Fig2, width = 10, height = 7)
ggsave(file = "Fig3.svg", plot = Fig3, width = 10, height = 24.5)
ggsave(file = "Fig4.svg", plot = Fig4, width = 10, height = 10.5)

ggsave(file = "Fig1.png", plot = Fig1, width = 10, height = 17.5)
ggsave(file = "Fig2.png", plot = Fig2, width = 10, height = 7)
ggsave(file = "Fig3.png", plot = Fig3, width = 10, height = 24.5)
ggsave(file = "Fig4.png", plot = Fig4, width = 10, height = 10.5)